
#include "femSolids.h"
#include "headersVTK.h"
#include "MyTime.h"
#include "log.h"
#include <string>
#include <boost/algorithm/string.hpp>
#include <sys/select.h>

extern  MyTime myTime;
extern  bool debug;

/// int ParMETIS_V3_PartMeshKway (
///     idx_t *elmdist, idx_t *eptr, idx_t *eind, idx_t *elmwgt, idx_t *wgtflag, idx_t *numflag,
///     idx_t *ncon, idx_t *ncommonnodes, idx_t *nparts, real t *tpwgts, real t *ubvec,
///     idx_t *options, idx_t *edgecut, idx_t *part, MPI Comm *comm
/// )
[[nodiscard]] PetscErrorCode femSolids::splitSolidElements() {

    LOG_INFO("femSolids::prepareInputData()  .... STARTED ... " 
        + string(__FILE__) + ":" + to_string(__LINE__));
    assert(ndim > 0 && ndim < 4);

    ndof = ndim;

    int idd = ElementTypeDataList[0]->getElemTypeNameNum();
    LOG_INFO("idd = " + to_string(idd));
    if( idd == ELEM_TRUSS_2D_NODES2 )
      ndof = 2;
    else if( idd == ELEM_TRUSS_3D_NODES2 )
      ndof = 3;
    if( (idd == ELEM_BEAM_2D_NODES2) || (idd == ELEM_BEAM_2D_NODES3) )
      ndof = 3;
    if( (idd == ELEM_BEAM_3D_NODES2) || (idd == ELEM_BEAM_3D_NODES3) )
      ndof = 6;
    else if( idd == ELEM_SHELL_FLAT_QUAD4 )
      ndof = 6;

    /*** parameter elmdist ***/
    LOG_INFO("prepare parameter elmdist " + string(__FILE__) + ":" + to_string(__LINE__));
    idx_t *elmdist;
    errpetsc = PetscMalloc1(n_mpi_procs + 1, &elmdist);CHKERRQ(errpetsc);
    LOG_INFO("PetscMalloc1 elmdist end " + string(__FILE__) + ":" + to_string(__LINE__));
    for (int i = 0; i < n_mpi_procs; i++) {
        elmdist[i] = i * ele_local_size;
    }
    elmdist[n_mpi_procs] = nElem_global;

    /*** parameter eptr,eind ***/
    LOG_INFO("prepare parameter eptr,eind " + string(__FILE__) + ":" + to_string(__LINE__));
    int npElem = elemConn[0].size() - 3;
    LOG_INFO("size of npElem = " + to_string(npElem));
    idx_t *eptr, *eind;
    errpetsc = PetscMalloc1(nElem_local + 1, &eptr);CHKERRQ(errpetsc);
    LOG_INFO("PetscMalloc1 eptr end "+ string(__FILE__) + ":" + to_string(__LINE__));
    errpetsc = PetscMalloc1(nElem_local * npElem, &eind);CHKERRQ(errpetsc);
    LOG_INFO("PetscMalloc1 eind end " + string(__FILE__) + ":" + to_string(__LINE__));

    vector<idx_t> vecTemp;

    eptr[0] = 0;
    int kk = 0;
    for (int ee = 0; ee < nElem_local; ee++) {
        eptr[ee + 1] = (ee + 1) * npElem;
        vecTemp = elemConn[ee];

        for(int ii = 0; ii < npElem; ii++) {
            eind[kk + ii] = vecTemp[3 + ii];
        }
        kk += npElem;
    }

    elemConn.clear();

    int ncommon_nodes;
    if(ndim == 2)
    {
      if((npElem == 3) || (npElem == 4))   // 3-noded tria or 4-noded quad elements
        ncommon_nodes = 2;
      else
        ncommon_nodes = 3;                 // 6-noded tria or 9-noded quad elements
    }
    else
    {
      if(npElem == 4)                      // 4-noded tetra element
        if(ndof > 3)
          ncommon_nodes = 2;
        else
          ncommon_nodes = 3;
      else if(npElem == 10)                // 10-noded tetra element
        ncommon_nodes = 6;
      else if(npElem == 6)                 // 6-noded Wedge/Prism element
        ncommon_nodes = 4;
      else if(npElem == 8)                 // 8-noded hexa element
        ncommon_nodes = 4;
      else if(npElem == 18)                // 18-noded Wedge/Prism element
        ncommon_nodes = 9;
      else
        ncommon_nodes = 9;                 // 27-noded hexa element
    }
    LOG_INFO("size of ncommon_nodes = " + to_string(ncommon_nodes) 
        + string(__FILE__) + ":" + to_string(__LINE__));

    /*** parameter others ***/
    LOG_INFO("prepare parameter others " + string(__FILE__) + ":" + to_string(__LINE__));
    idx_t *elmwgt = NULL;

    idx_t wgtflag_[] = {0};
    idx_t *wgtflag = wgtflag_;

    idx_t numflag_[] = {0};
    idx_t *numflag = numflag_;

    idx_t ncon_[] = {1};
    idx_t *ncon = ncon_;

    idx_t ncommonnodes_[] = {ncommon_nodes};
    idx_t *ncommonnodes = ncommonnodes_;

    idx_t nparts_[] = {n_mpi_procs};
    idx_t *nparts = nparts_;

    real_t *tpwgts;
    errpetsc = PetscMalloc1(n_mpi_procs * ncon[0], &tpwgts);CHKERRQ(errpetsc);
    LOG_INFO("PetscMalloc1 tpwgts end "+ string(__FILE__) + ":" + to_string(__LINE__));
    for (int i = 0; i < n_mpi_procs * ncon[0]; i++) {
        tpwgts[i] = 1.0f / n_mpi_procs;
    }

    real_t ubvec_[] = {1.05};
    real_t *ubvec = ubvec_;

    idx_t options_[METIS_NOPTIONS];
    METIS_SetDefaultOptions(options_);
    options_[METIS_OPTION_PTYPE]     = METIS_PTYPE_KWAY;
    options_[METIS_OPTION_OBJTYPE]   = METIS_OBJTYPE_VOL;
    options_[METIS_OPTION_NUMBERING] = 0;
    idx_t *options = options_;

    idx_t *edgecut;
    PetscMalloc1(1, &edgecut);CHKERRQ(errpetsc);
    LOG_INFO("PetscMalloc1 edgecut end "+ string(__FILE__) + ":" + to_string(__LINE__));

    // idx_t *part
    PetscMalloc1(nElem_local, &part);CHKERRQ(errpetsc);
    LOG_INFO("PetscMalloc1 part end "+ string(__FILE__) + ":" + to_string(__LINE__));

    MPI_Comm comm = static_cast<MPI_Comm>(MPI_COMM_WORLD);

    /*** prepare data end ***/
    LOG_INFO("prepare data end "+ string(__FILE__) + ":" + to_string(__LINE__));

    int parmetis_ret = ParMETIS_V3_PartMeshKway(
        elmdist, eptr, eind, elmwgt, wgtflag, numflag, ncon, ncommonnodes,
        nparts, tpwgts, ubvec, options, edgecut, part, &comm
    );

    if(parmetis_ret == METIS_OK) {
        LOG_INFO("PARMETIS partition routine successful " 
            + string(__FILE__) + ":" + to_string(__LINE__));
    } else {
        LOG_ERROR("PARMETIS partition routine FAILED, err_code: " + to_string(parmetis_ret)
            + string(__FILE__) + ":" + to_string(__LINE__));
        exit(1);
    }
    
    errpetsc = PetscFree(eptr  ); CHKERRQ(errpetsc);
    errpetsc = PetscFree(eind  ); CHKERRQ(errpetsc);
    errpetsc = PetscFree(tpwgts); CHKERRQ(errpetsc);
    
    MPI_Barrier(MPI_COMM_WORLD);

    for (int i = 0; i < nElem_local; i++) {
        LOG_DEBUG("ele[" + to_string(nElem_start + i) + "]\t --> " + to_string(part[i]));
    }

    return errpetsc;
}

[[nodiscard]] PetscErrorCode femSolids::recvLocalElements() {

    vector<int> pid_ele_count(FEM_MAX_MPI_NUM, 0);
    vector<int> pid_ele_idx[FEM_MAX_MPI_NUM];

    for (int i = 0; i < nElem_local; i++) {
        int target_pid = part[i];
        int ele_global_idx = nElem_start + i;
        pid_ele_count[target_pid]++;
        pid_ele_idx[target_pid].push_back(ele_global_idx);
    }

    vector<int> part_count(FEM_MAX_MPI_NUM, 0);
    vector<int> part_all[FEM_MAX_MPI_NUM];
    if(this_mpi_proc == 0) {
        ele_idx_map.resize(nElem_global);
    }
  
    if(this_mpi_proc != 0) {
        MPI_Send(&nElem_local, 1, MPI_INT, 0, MPI_ELE_JOIN_NUM, MPI_COMM_WORLD);
    }

    if(this_mpi_proc == 0) {
        for(int pid_send = 1; pid_send < n_mpi_procs; pid_send++) {
            MPI_Recv(&part_count[pid_send], 1, MPI_INT, pid_send, MPI_ELE_JOIN_NUM, MPI_COMM_WORLD, MPI_STATUS_IGNORE);
            part_all[pid_send].resize(part_count[pid_send]);
            LOG_INFO("recv part size: " + to_string(part_count[pid_send]) + " from pid: " + to_string(pid_send));
        }
    }
    MPI_Barrier(MPI_COMM_WORLD);
  
    if(this_mpi_proc != 0) {
        MPI_Send(&part[0], nElem_local, MPI_INT, 0, MPI_ELE_JOIN, MPI_COMM_WORLD);
    }

    if(this_mpi_proc == 0) {
        int base_idx;
        for(int pid_send = 1; pid_send < n_mpi_procs; pid_send++) {
            LOG_INFO("recv part start from pid: " + to_string(pid_send));
            MPI_Recv(part_all[pid_send].data(), part_count[pid_send], MPI_INT, pid_send, MPI_ELE_JOIN, MPI_COMM_WORLD, MPI_STATUS_IGNORE);
            base_idx = pid_send * ele_local_size;
            for(int i = 0; i < part_count[pid_send]; i++) {
                ele_idx_map[base_idx + i] = part_all[pid_send][i];
            }
            LOG_INFO("recv part end from pid: " + to_string(pid_send));
        }
        base_idx = 0;
        for(int i = 0; i < nElem_local; i++) {
           ele_idx_map[base_idx + i] = part[i];
        }

        for(int i = 0; i < nElem_global; i++) {
            LOG_DEBUG("ele[" + to_string(i) + "]\t-->" + to_string(ele_idx_map[i]));
        }
    }
    LOG_INFO("recv part end: " + string(__FILE__) + ":" + to_string(__LINE__));
    MPI_Barrier(MPI_COMM_WORLD);

    errpetsc = PetscFree(part); CHKERRQ(errpetsc);

    /*** Send ele count ***/
    for (int pid_recv = 0; pid_recv < n_mpi_procs; pid_recv++) {
        if(pid_recv != this_mpi_proc) {
            MPI_Send(&pid_ele_count[pid_recv], 1, MPI_INT, pid_recv, MPI_SEND_ELE_NUM, MPI_COMM_WORLD);

            LOG_INFO("Send ele from pid[" + to_string(this_mpi_proc) + "] to pid[" + to_string(pid_recv) +
                "]\tcount: " + to_string(pid_ele_count[pid_recv])
                + string(__FILE__) + ":" + to_string(__LINE__));
        }
    }

    /*** Recv ele count ***/
    vector<int> pid_ele_recv(FEM_MAX_MPI_NUM, 0);

    for (int pid_send = 0; pid_send < n_mpi_procs; pid_send++) {
        if(pid_send != this_mpi_proc) {
            MPI_Recv(&pid_ele_recv[pid_send], 1, MPI_INT, pid_send, MPI_RECV_ELE_NUM, MPI_COMM_WORLD, MPI_STATUS_IGNORE);

            LOG_INFO("Recv ele from pid[" + to_string(pid_send) + "] to pid[" + to_string(this_mpi_proc) +
                "]\tcount: " + to_string(pid_ele_recv[pid_send])
                + string(__FILE__) + ":" + to_string(__LINE__));
        }
    }

    /*** Init recv array ***/
    vector<int> pid_ele_idx_recv[FEM_MAX_MPI_NUM];

    for (int pid_send = 0; pid_send < n_mpi_procs; pid_send++) {
        if(pid_send != this_mpi_proc) {
            pid_ele_idx_recv[pid_send].resize(pid_ele_recv[pid_send]);
            LOG_INFO("Recv ele idx_array[" + to_string(pid_send) + "]\tresize: " + to_string(pid_ele_idx_recv[pid_send].size())
                + string(__FILE__) + ":" + to_string(__LINE__));
        }
    }
    LOG_INFO("Pid_ele_idx_recv init end " + string(__FILE__) + ":" + to_string(__LINE__));

    MPI_Barrier(MPI_COMM_WORLD);

    vector<MPI_Request> send_requests(n_mpi_procs);
    vector<MPI_Request> recv_requests(n_mpi_procs);

    /*** Send ele start ***/
    for (int pid_recv = 0; pid_recv < n_mpi_procs; pid_recv++) {
       if(pid_recv != this_mpi_proc) {
           //MPI_Send(pid_ele_idx[pid_recv].data(), pid_ele_count[pid_recv], MPI_INT, pid_recv, MPI_SEND_ELE, MPI_COMM_WORLD);
           MPI_Isend(pid_ele_idx[pid_recv].data(), pid_ele_count[pid_recv], MPI_INT, pid_recv, MPI_SEND_ELE, MPI_COMM_WORLD, &send_requests[pid_recv]);

           LOG_INFO("Send ele start, form pid[" + to_string(this_mpi_proc) + "] to pid[" + to_string(pid_recv) + "]"
               + string(__FILE__) + ":" + to_string(__LINE__));
           LOG_INFO("Send size: " + to_string(pid_ele_count[pid_recv]) + 
               ", Array size: " + to_string(pid_ele_idx[pid_recv].size()));
       }
    }

    /*** Recv ele start ***/
    for (int pid_send = 0; pid_send < n_mpi_procs; pid_send++) {
        if(pid_send != this_mpi_proc) {
            //MPI_Recv(pid_ele_idx_recv[pid_send].data(), pid_ele_recv[pid_send], MPI_INT, pid_send, MPI_RECV_ELE, MPI_COMM_WORLD, MPI_STATUS_IGNORE);
            MPI_Irecv(pid_ele_idx_recv[pid_send].data(), pid_ele_recv[pid_send], MPI_INT, pid_send, MPI_RECV_ELE, MPI_COMM_WORLD, &recv_requests[pid_send]);

            LOG_INFO("Recv ele start, from pid[" + to_string(pid_send) + "] to pid[" + to_string(this_mpi_proc) + "]"
                + string(__FILE__) + ":" + to_string(__LINE__));
            LOG_INFO("Recv size: " + to_string(pid_ele_recv[pid_send]) + 
               ", Array size: " + to_string(pid_ele_idx_recv[pid_send].size()));
        }
    }

    MPI_Waitall(n_mpi_procs, send_requests.data(), MPI_STATUSES_IGNORE);
    MPI_Waitall(n_mpi_procs, recv_requests.data(), MPI_STATUSES_IGNORE);

    for (int pid_send = 0; pid_send < n_mpi_procs; pid_send++) {
        if (pid_send != this_mpi_proc) {
            MPI_Wait(&recv_requests[pid_send], MPI_STATUS_IGNORE);
        }
    }

    LOG_INFO("Recv ele end " + string(__FILE__) + ":" + to_string(__LINE__));

    for (int pid = 0; pid < n_mpi_procs; pid++) {
        if(pid != this_mpi_proc) {
            for (int i = 0; i < pid_ele_recv[pid]; i++) {
                ele_local_idx_set.insert(pid_ele_idx_recv[pid][i]);
            }
        } else {
            for (int i = 0; i < pid_ele_count[pid]; i++) {
                ele_local_idx_set.insert(pid_ele_idx[pid][i]);
            }
        } 
    }

    MPI_Barrier(MPI_COMM_WORLD);
    LOG_INFO("Insert ele set end " + string(__FILE__) + ":" + to_string(__LINE__));

    for (auto &idx : ele_local_idx_set) {
        LOG_DEBUG("ele[" + to_string(idx) + "]\t--> " + to_string(this_mpi_proc));
    }

    return errpetsc;
}

int  femSolids::postProcess()
{
    if(this_mpi_proc != 0)
      return 0;

    if(debug) PetscPrintf(MPI_COMM_WORLD, " \n femSolids::postProcess \n");

    struct timeval start, end;
    gettimeofday(&start, NULL);
    //
    // setup and write vtk data
    //
    //////////////////////////////////////////////


    vtkSmartPointer<vtkUnstructuredGrid>     uGridVTK     =  vtkSmartPointer<vtkUnstructuredGrid>::New();
    vtkSmartPointer<vtkPoints>               pointsVTK    =  vtkSmartPointer<vtkPoints>::New();
    vtkSmartPointer<vtkVertex>               vertexVTK    =  vtkSmartPointer<vtkVertex>::New();

    vtkSmartPointer<vtkLine>                 lineVTK      =  vtkSmartPointer<vtkLine>::New();
    vtkSmartPointer<vtkTriangle>             triaVTK      =  vtkSmartPointer<vtkTriangle>::New();
    vtkSmartPointer<vtkQuad>                 quadVTK      =  vtkSmartPointer<vtkQuad>::New();
    vtkSmartPointer<vtkTetra>                tetraVTK     =  vtkSmartPointer<vtkTetra>::New();
    vtkSmartPointer<vtkWedge>                wedgeVTK     =  vtkSmartPointer<vtkWedge>::New();
    vtkSmartPointer<vtkHexahedron>           hexaVTK      =  vtkSmartPointer<vtkHexahedron>::New();

    vtkSmartPointer<vtkQuadraticTriangle>    tria6VTK     =  vtkSmartPointer<vtkQuadraticTriangle>::New();
    vtkSmartPointer<vtkBiQuadraticQuad>      quad9VTK     =  vtkSmartPointer<vtkBiQuadraticQuad>::New();
    vtkSmartPointer<vtkQuadraticTetra>       tetra10VTK   =  vtkSmartPointer<vtkQuadraticTetra>::New();
    vtkSmartPointer<vtkBiQuadraticQuadraticWedge>  wedge18VTK = vtkSmartPointer<vtkBiQuadraticQuadraticWedge>::New();
    vtkSmartPointer<vtkQuadraticHexahedron>  hexa27VTK    =  vtkSmartPointer<vtkQuadraticHexahedron>::New();

    vtkSmartPointer<vtkFloatArray>           vecVTK       =  vtkSmartPointer<vtkFloatArray>::New();
    vtkSmartPointer<vtkFloatArray>           vecVTK2      =  vtkSmartPointer<vtkFloatArray>::New();
    vtkSmartPointer<vtkFloatArray>           scaVTK       =  vtkSmartPointer<vtkFloatArray>::New();
    vtkSmartPointer<vtkFloatArray>           scaVTK2      =  vtkSmartPointer<vtkFloatArray>::New();

    vtkSmartPointer<vtkFloatArray>           cellDataVTK  =  vtkSmartPointer<vtkFloatArray>::New();
    vtkSmartPointer<vtkFloatArray>           cellDataVTK2 =  vtkSmartPointer<vtkFloatArray>::New();
    vtkSmartPointer<vtkFloatArray>           strainVTK    =  vtkSmartPointer<vtkFloatArray>::New();
    vtkSmartPointer<vtkFloatArray>           stressVTK    =  vtkSmartPointer<vtkFloatArray>::New();
    vtkSmartPointer<vtkFloatArray>           procIdVTKcell =  vtkSmartPointer<vtkFloatArray>::New();

    vtkSmartPointer<vtkFloatArray>           timeloadstamp  =  vtkSmartPointer<vtkFloatArray>::New();

    vtkSmartPointer<vtkXMLUnstructuredGridWriter>  writerUGridVTK =  vtkSmartPointer<vtkXMLUnstructuredGridWriter>::New();


    int  dd, ii, jj, kk, ll, nn, nlocal, ind1, ind2, e, ee, count, gcount, ind, npElem;
    int  n1, n2, n3, n4, n5, n6;
    double vec[3]={0.0,0.0,0.0}, vec2[3]={0.0,0.0,0.0}, xx, yy, zz, val, maxdisp = 0.0, maxstress = 0.0;
    //                        {0,1,2,3,4,5,6,7,8,9, 10,11,12,13,14,15,16,17,18,19,20,21,22,23,24,25,26};
    int  hexa27nodemap[27] =  {0,1,2,3,4,5,6,7,8,11,16, 9,17,10,18,19,12,15,13,14,24,22,20,21,23,25,26};
    //                        {0,1,2,3,4,5,6,7, 8,9,10,11,12,13,14,15,16,17};
    int  wedge18nodemap[18] = {0,1,2,3,4,5,6,8,12,7,13,14, 9,11,10,15,17,16};
    //                        {0,1,2,3,4,5,6,7,8,9};
    int  tetra10nodemap[10] = {0,1,2,3,4,5,6,7,9,8};

    int vartype=0, varindex=0, index=0;
    bool extrapolateFlag=true;
    vector<int>  nodeNums;

    vtkIdType pt[50];

    vecVTK->SetName("disp");
    //vecVTK2->SetName("velo");
    //if(ndof > 1)
    //{
      vecVTK->SetNumberOfComponents(3);
      vecVTK->SetNumberOfTuples(nNode_global);
      //vecVTK2->SetNumberOfComponents(3);
      //vecVTK2->SetNumberOfTuples(nNode_global);
    //}

    procIdVTKcell->SetName("procID");
    procIdVTKcell->SetNumberOfTuples(nElem_global);

    timeloadstamp->SetName("TIME");
    timeloadstamp->SetNumberOfTuples(1);
    if(ARC_LENGTH)
      timeloadstamp->SetTuple1(0, loadFactor);
    else
      timeloadstamp->SetTuple1(0, myTime.cur);


    //int idd = SolnData.ElemProp[0].id;
    int idd = 2;

    ifstream infile;
    infile.open(input_file_name);
    if(infile.fail()) {
        LOG_ERROR("Data file open failed" + string(__FILE__) + ":" + to_string(__LINE__));
        exit(1);
    }

    LOG_INFO("join output start " + string(__FILE__) + ":" + to_string(__LINE__));

    if(ndim == 2)
    {
      for(ii=0;ii<nNode_global;ii++)
      {
          nn = node_map_get_new[ii];

          xx = GeomData.NodePosNew[nn][0];
          yy = GeomData.NodePosNew[nn][1];
          
          pt[0] = pointsVTK->InsertNextPoint(xx, yy, 0.0);

          vertexVTK->GetPointIds()->SetId(0, pt[0]);

          kk = nn*ndof;

          vec[0] = SolnData.disp[kk];
          vec[1] = SolnData.disp[kk+1];

          vecVTK->InsertTuple(ii, vec);

      }

      string line;

      while(getline(infile,line)) {
          boost::trim(line);
          if(isAValidLine(line)) {
              if((line.compare(string("Elements")) == 0) || (line.compare(string("elements")) == 0)) {
                  getline(infile,line);    

                  vector<string>  stringlist;
                  vector<int>     nodeNums;
                  for(ee = 0; ee < nElem_global; ee++) {
                      getline(infile,line);
                      boost::trim(line);
                      boost::algorithm::split(stringlist, line, boost::is_any_of("\t "), boost::token_compress_on);

                      npElem = stringlist.size() - 4;
                      nodeNums.resize(npElem);
                      for(int i = 0; i < npElem; i++) {
                         nodeNums[i] = stoi(stringlist[i + 4]) - 1;
                      }

                      procIdVTKcell->SetTuple1(ee, ele_idx_map[ee]);
                      LOG_DEBUG("join[" + to_string(ee) + "]\t--> " + to_string(ele_idx_map[ee]));

                      if(npElem == 2)
                      {
                        lineVTK->GetPointIds()->SetId(0, node_map_get_old[nodeNums[0]]);
                        lineVTK->GetPointIds()->SetId(1, node_map_get_old[nodeNums[1]]);

                        uGridVTK->InsertNextCell(lineVTK->GetCellType(), lineVTK->GetPointIds());
                      }
                      else if(npElem == 3)
                      {
                        for(ii=0; ii<npElem; ii++)
                          triaVTK->GetPointIds()->SetId(ii, node_map_get_old[nodeNums[ii]]);
                        
                        uGridVTK->InsertNextCell(triaVTK->GetCellType(), triaVTK->GetPointIds());
                      }
                      else if(npElem == 6)
                      {
                        for(ii=0; ii<npElem; ii++)
                          tria6VTK->GetPointIds()->SetId(ii, nodeNums[ii]);
                        
                        uGridVTK->InsertNextCell(tria6VTK->GetCellType(), tria6VTK->GetPointIds());
                      }
                      else if(npElem == 4)
                      {
                        for(ii=0; ii<npElem; ii++)
                          quadVTK->GetPointIds()->SetId(ii, node_map_get_old[nodeNums[ii]]);
                        
                        uGridVTK->InsertNextCell(quadVTK->GetCellType(), quadVTK->GetPointIds());
                      }
                      else if(npElem == 9)
                      {
                        for(ii=0; ii<npElem; ii++)
                          quad9VTK->GetPointIds()->SetId(ii, nodeNums[ii]);
                        
                        uGridVTK->InsertNextCell(quad9VTK->GetCellType(), quad9VTK->GetPointIds());
                      }
                      else
                      {
                        PetscPrintf(MPI_COMM_WORLD, "\n\n Wrong 2D element type for post-processing \n\n");
                        exit(1);
                      }
                  }
              }
          }
      }
    }
    else //if(ndim == 3)
    {
      for(ii=0;ii<nNode_global;ii++)
      {
          nn = node_map_get_new[ii];

          xx = GeomData.NodePosNew[nn][0];
          yy = GeomData.NodePosNew[nn][1];
          zz = GeomData.NodePosNew[nn][2];

          pt[0] = pointsVTK->InsertNextPoint(xx, yy, zz);

          vertexVTK->GetPointIds()->SetId(0, pt[0]);

          kk = nn*ndof;

          vec[0] = SolnData.disp[kk];
          vec[1] = SolnData.disp[kk+1];
          vec[2] = SolnData.disp[kk+2];

          printf("vec0: %.20f vec1: %.20f vec2: %.20f\n",
              vec[0], vec[1], vec[2]);

          vecVTK->InsertTuple(ii, vec);

      }
      string line;

      while(getline(infile,line)) {
          boost::trim(line);
          if(isAValidLine(line)) {
              if((line.compare(string("Elements")) == 0) || (line.compare(string("elements")) == 0)) {
                  getline(infile,line);    

                  vector<string>  stringlist;
                  vector<int>     nodeNums;
                  for(ee = 0; ee < nElem_global; ee++) {
                      getline(infile,line);
                      boost::trim(line);
                      boost::algorithm::split(stringlist, line, boost::is_any_of("\t "), boost::token_compress_on);

                      npElem = stringlist.size() - 4;
                      nodeNums.resize(npElem);
                      for(int i = 0; i < npElem; i++) {
                         nodeNums[i] = stoi(stringlist[i + 4]) - 1;
                      }

                      procIdVTKcell->SetTuple1(ee, ele_idx_map[ee]);
                      LOG_DEBUG("join[" + to_string(ee) + "]\t--> " + to_string(ele_idx_map[ee]));

                      if(npElem == 4) // tet
                      {
                        if(elems[0]->ElemTypeData->getElemTypeNameNum() == ELEM_SOLID_3D_TETR4)
                        {
                          for(ii=0; ii<npElem; ii++)
                            tetraVTK->GetPointIds()->SetId(ii, nodeNums[ii]);
                          
                          uGridVTK->InsertNextCell(tetraVTK->GetCellType(), tetraVTK->GetPointIds());
                        }
                        else
                        {
                          for(ii=0; ii<npElem; ii++)
                            quadVTK->GetPointIds()->SetId(ii, nodeNums[ii]);
                          
                          uGridVTK->InsertNextCell(quadVTK->GetCellType(), quadVTK->GetPointIds());
                        }
                      }
                      else if(npElem == 6)
                      {
                        for(ii=0; ii<npElem; ii++)
                          wedgeVTK->GetPointIds()->SetId(ii, nodeNums[ii]);
                        
                        uGridVTK->InsertNextCell(wedgeVTK->GetCellType(), wedgeVTK->GetPointIds());
                      }
                      else if(npElem == 8) // hexa
                      {
                        for(ii=0; ii<npElem; ii++)
                          hexaVTK->GetPointIds()->SetId(ii, nodeNums[ii]);
                        
                        uGridVTK->InsertNextCell(hexaVTK->GetCellType(), hexaVTK->GetPointIds());
                      }
                      else if(npElem == 10) // 10-node tetrahedron
                      {
                        for(ii=0; ii<npElem; ii++)
                          tetra10VTK->GetPointIds()->SetId(tetra10nodemap[ii], nodeNums[ii]);
                        
                        uGridVTK->InsertNextCell(tetra10VTK->GetCellType(), tetra10VTK->GetPointIds());
                      }
                      else if(npElem == 18) // 18-node wedge
                      {
                        for(ii=0; ii<npElem; ii++)
                          wedge18VTK->GetPointIds()->SetId(wedge18nodemap[ii], nodeNums[ii]);
                        
                        uGridVTK->InsertNextCell(wedge18VTK->GetCellType(), wedge18VTK->GetPointIds());
                      }
                      else if(npElem == 27) // 27-node hexahedron
                      {
                        for(ii=0; ii<npElem; ii++)
                          wedge18VTK->GetPointIds()->SetId(hexa27nodemap[ii], nodeNums[ii]);
                        
                        uGridVTK->InsertNextCell(wedge18VTK->GetCellType(), wedge18VTK->GetPointIds());
                      }
                      else
                      {
                        PetscPrintf(MPI_COMM_WORLD, "\n\n Wrong 3D element type for post-processing \n\n");
                        exit(1);
                      }
                  }
              }
          }
      }
    }

    LOG_INFO("join output end " + string(__FILE__) + ":" + to_string(__LINE__));
    infile.close();

    //cout << " vartype " << endl;
    //cout << vartype << '\t' << varindex << '\t' << index << endl;

    cellDataVTK->SetName("pres");
    cellDataVTK->SetNumberOfTuples(nElem_global);
    for(ee=0;ee<nElem_global;ee++)
    {
      //TODO
      //elems[ee]->elementContourplot(vartype, varindex, index);

      //cellDataVTK->InsertTuple1(ee, elems[ee]->vals2project[0]);

      //maxstress = max(maxstress, elems[ee]->vals2project[0]);
    }

    if(intVarFlag)
    {
      scaVTK2->SetName("epstrn");
      scaVTK2->SetNumberOfTuples(count);
      projectStresses(1, 5, 6, index);
    }

    scaVTK->SetName("pres");
    scaVTK->SetNumberOfTuples(count);

    //cout << vartype << '\t' << varindex << '\t' << index << endl;

    //projectStresses(extrapolateFlag, vartype, varindex, index);
    //cout << vartype << '\t' << varindex << '\t' << index << endl;
    //projectFromElemsToNodes(extrapolateFlag, vartype, varindex, index);


    uGridVTK->SetPoints(pointsVTK);
    uGridVTK->GetFieldData()->AddArray(timeloadstamp);
    //uGridVTK->GetPointData()->AddArray(scaVTK);
    uGridVTK->GetPointData()->AddArray(vecVTK);
    //uGridVTK->GetPointData()->AddArray(vecVTK2);
    uGridVTK->GetCellData()->AddArray(procIdVTKcell);
    //if(intVarFlag)
      //uGridVTK->GetPointData()->AddArray(scaVTK2);


    //Write the file
    char VTKfilename[200];
    sprintf(VTKfilename,"%s%s%06d%s", jobname.c_str(), "-",fileCount,".vtu");

    writerUGridVTK->SetFileName(VTKfilename);
    writerUGridVTK->SetInputData(uGridVTK);
    writerUGridVTK->Write();

    gettimeofday(&end, NULL);
    double write_time = get_time_fem(&start, &end);
    LOG_INFO("write data " + to_string(fileCount) + " end cost time(s): " + to_string(write_time));

    fileCount++;

    return 0;
}



int  femSolids::projectStrains(bool extrapolateFlag, int vartype, int varindex, int index)
{
    return 0;
}



int  femSolids::projectStresses(bool extrapolateFlag, int vartype, int varindex, int index)
{
    int  ee, ii;

    //cout << " extrapolateFlag = " << extrapolateFlag << endl;

    // compute the stresses at element nodes
    for(ee=0; ee<nElem_global; ee++)
    {
      elems[ee]->projectToNodes(extrapolateFlag, vartype, varindex, index);
    }

    vector<int>  nodeNums;
    vector<double>  output(nNode_global, 0.0);
    double  volu, val;

    for(ee=0; ee<nElem_global; ee++)
    {
      nodeNums = elems[ee]->nodeNums;

      elems[ee]->computeVolume(false);

      volu = elems[ee]->getVolume();

      //cout << " volu= " << volu << endl;

      for(ii=0; ii<nodeNums.size(); ii++)
      {
        //output[nodeNums[ii]] += volu*elems[ee]->vals2project[ii];
        output[nodeNums[ii]] += elems[ee]->vals2project[ii];
      }
    }

    for(ee=0; ee<nNode_global; ee++)
    {
      //printVector(node_elem_conn[ee]);
      volu = 0.0;
      for(ii=0; ii<node_elem_conn[ee].size(); ii++)
      {
        volu += elems[node_elem_conn[ee][ii]]->getVolume();
      }

      //output[ee] /= volu;
      output[ee] /= node_elem_conn[ee].size();
    }

    // for Bernstein elements, the values for midnodes need to be interpolated
    // Only to be done when the values are extrapolated from Gauss points
    //if(extrapolateFlag)
    //{
      for(ii=0; ii<nNode_global; ii++)
      {
        if( midnodeData[ii][0] )
        {
          val  = 0.50*output[ii];
          val += 0.25*output[midnodeData[ii][1]];
          val += 0.25*output[midnodeData[ii][2]];

          //scaVTK->SetTuple1(ii, val);
        }
        //else
          //scaVTK->SetTuple1(ii, output[ii]);
      }
    //}

    //cout << " extrapolateFlag = " << extrapolateFlag << endl;

    return 0;
}


int  femSolids::projectInternalVariables(bool extrapolateFlag, int vartype, int varindex, int index)
{
    return 0;
}






/*
int  femSolids::postProcess()
{
    cout << " femSolids::postProcess " << endl;

    //
    // setup and write vtk data
    //
    //////////////////////////////////////////////


    vtkSmartPointer<vtkUnstructuredGrid>     uGridVTK     =  vtkSmartPointer<vtkUnstructuredGrid>::New();
    vtkSmartPointer<vtkPoints>               pointsVTK    =  vtkSmartPointer<vtkPoints>::New();
    vtkSmartPointer<vtkVertex>               vertexVTK    =  vtkSmartPointer<vtkVertex>::New();

    vtkSmartPointer<vtkLine>                 lineVTK      =  vtkSmartPointer<vtkLine>::New();

    vtkSmartPointer<vtkTriangle>             triaVTK      =  vtkSmartPointer<vtkTriangle>::New();
    vtkSmartPointer<vtkQuad>                 quadVTK      =  vtkSmartPointer<vtkQuad>::New();
    vtkSmartPointer<vtkTetra>                tetraVTK     =  vtkSmartPointer<vtkTetra>::New();
    vtkSmartPointer<vtkHexahedron>           hexaVTK      =  vtkSmartPointer<vtkHexahedron>::New();

    vtkSmartPointer<vtkQuadraticTriangle>    tria6VTK     =  vtkSmartPointer<vtkQuadraticTriangle>::New();
    vtkSmartPointer<vtkBiQuadraticQuad>      quad9VTK     =  vtkSmartPointer<vtkBiQuadraticQuad>::New();
    vtkSmartPointer<vtkQuadraticTetra>       tetra10VTK   =  vtkSmartPointer<vtkQuadraticTetra>::New();
    //vtkSmartPointer<vtkQuadraticHexahedron>  hexa27VTK    =  vtkSmartPointer<vtkQuadraticHexahedron>::New();

    vtkSmartPointer<vtkFloatArray>           vecVTK       =  vtkSmartPointer<vtkFloatArray>::New();
    vtkSmartPointer<vtkFloatArray>           scaVTK       =  vtkSmartPointer<vtkFloatArray>::New();

    vtkSmartPointer<vtkXMLUnstructuredGridWriter>  writerUGridVTK =  vtkSmartPointer<vtkXMLUnstructuredGridWriter>::New();


    int  ee, ii, jj, kk, nn, n1, n2, npElem;
    double  val, vec[3]={0.0,0.0,0.0}, dsp[3]={0.0,0.0,0.0};

    vecVTK->SetNumberOfTuples(nNode_global_global);
    vecVTK->SetNumberOfComponents(3);
    scaVTK->SetNumberOfTuples(nNode_global_global);

    vecVTK->SetName("Disp");

    vtkIdType pt[10];

    if(ndim == 2)
    {
      for(ii=0; ii<nNode_global_global; ii++)
      {
        pt[0] = pointsVTK->InsertNextPoint(nodeCoordsCur[ii][0], nodeCoordsCur[ii][1], 0.0);

        kk = ii*ndof;

        dsp[0] = SolnData.disp[kk];
        dsp[1] = SolnData.disp[kk+1];

        vecVTK->InsertTuple(ii, dsp);
      }

      for(ee=0; ee<nElem_global_global; ee++)
      {
        npElem = elems[ee]->nodeNums.size();

        if( npElem == 2 )
        {
          for(ii=0; ii<npElem; ii++)
            lineVTK->GetPointIds()->SetId(ii, elems[ee]->nodeNums[ii] );

          uGridVTK->InsertNextCell(lineVTK->GetCellType(), lineVTK->GetPointIds());
        }
        else if(npElem == 3)
        {
          vecVTK->InsertTuple(elemConn[ee][8], vec);

          uGridVTK->InsertNextCell(quad9VTK->GetCellType(), quad9VTK->GetPointIds());
        }
      }
    }
    else // (ndim == 3)
    {
      for(ii=0; ii<nNode_global_global; ii++)
      {
        pt[0] = pointsVTK->InsertNextPoint(nodeCoordsOrig[ii][0], nodeCoordsOrig[ii][1], nodeCoordsOrig[ii][2]);
      }

      for(ee=0; ee<nElem_global_global; ee++)
      {
        npElem = elemConn[ee].size();
      }
    }

    uGridVTK->SetPoints(pointsVTK);
    //uGridVTK->GetPointData()->SetScalars(scaVTK);
    uGridVTK->GetPointData()->SetVectors(vecVTK);

    //Write the file.

    char VTKfilename[200];

    sprintf(VTKfilename,"%s%s%06d%s", dirname.c_str(), "-solid-",fileCount,".vtu");

    writerUGridVTK->SetFileName(VTKfilename);

    writerUGridVTK->SetInputData(uGridVTK);
    writerUGridVTK->Write();

    fileCount = fileCount+1;

    return 0;
}
*/
